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ABSTRACT 

Recent simulations of self-gravitating accretion discs, carried out using a three-dimensional 
Smoothed Particle Hydrodynamics (SPH) code by Meru and Bate, have been interpreted as 
implying that three-dimensional global discs fragment much more easily than would be ex- 
pected from a two-dimensional local model. Subsequently, global and local two-dimensional 
models have been shown to display similar fragmentation properties, leaving it unclear 
whether the three-dimensional results reflect a physical effect or a numerical problem as- 
sociated with the treatment of cooling or artificial viscosity in SPH. Here, we study how frag- 
mentation of self-gravitating disc flows in SPH depends upon the implementation of cooling. 
We run disc simulations that compare a simple cooling scheme, in which each particle loses 
energy based upon its internal energy per unit mass, with a method in which the cooling is 
derived from a smoothed internal energy density field. For the simple per particle cooling 
scheme, we find a significant increase in the minimum cooling time scale for fragmentation 
with increasing resolution, matching previous results. Switching to smoothed cooling, how- 
ever, results in lower critical cooling time scales, and tentative evidence for convergence at 
the highest spatial resolution tested. We conclude that precision studies of fragmentation using 
SPH require careful consideration of how cooling (and, probably, artificial viscosity) is imple- 
mented, and that the apparent non-convergence of the fragmentation boundary seen in prior 
simulations is likely a numerical effect. In real discs, where cooling is physically smoothed 
by radiative transfer effects, the fragmentation boundary is probably displaced from the two- 
dimensional value by a factor that is only of the order of unity. 
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1 INTRODUCTION 

Accretion discs exist in many astrophysical systems - in x-ray bi- 
naries and cataclysmic variable stars, around supermassive black 
holes in the nuclei of active galaxies, and around young, newly 
forming stars. In some cases, in particular in active galactic nu- 
clei and around young stars, these discs can be sufficiently massive 
that their own self-gravity plays an important role in their evolu- 
tion. The susceptibilty of an infinitesimally thin accretion disc to 
the growth of the fflavitat ional instabilit y can be determined using 
the Toomre Q parameter (iToomrell 19641) 



where Cs is the disc sound speed, Q is the angular frequency, G is 
the gravitational constant, and S is the disc surface mass density. 



E-mail: wknir@roe.ac.uk 



Discs are unstable to a xisymmetric perturb ations if Q < 1, while 
numerical simulations ( Durisen et al.l2007 l) suggest that the growth 
of non-axisymmetric perturbations occurs for Q < 1.5 — 1.7. 

It has been understood for quite some time that the 
gravitational instability can act to transport angular momen- 
tumm outwards, allowing ma ss to accrete onto the cen- 
tral object jLin & Pringle 1987:1: iLaughlin & Bodenheimej Il994 
lArmitage. Livio & Pringle..200lh and could well be the primary 
transport mechanisni during the earliest stages of sta r formation 
iUn & Pringli 1 1981 iRice. Mavo & Armitage 1 12010). It is also 
possible that, if sufficie ntly unstable, these discs may fragment 
to form bound objects (lKuipei]ll95lh . In the case of protostel- 
lar discs, these ob jects could contract to form gas giant planets 
(IBossII 19981 . [200Qh or could form stars in discs around supermas- 
sive bl ack holes ( Sh losman & BegelmanI Il989l : iGoodmanI l2003l : 
iBonnell & Rice 2008,) . 

IPaczvnskil (Il978l) suggested that non-fragmenti ng, self- 
gravitating discs may exist in a state of marginal stability [Gammid 
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(l200lh . using two-dimensional, shearing sh eet simulations , quanti- 
fied this by showing that, in addition to the lXoomr^ (Il964h param- 
eter, the evolution of a self-gravitating accretion disc is also con- 
trolled by the rate at which it loses energy. Using a cooling time, 
Tc, of the following form 

Tc^P^~\ (2) 

where /3 is a constant, 'Gammie" ("2001) showed that if /3 > 3, the 
system settles into a quasi-steady, marginally stable state. If not, the 
system fragments to for m bou nd objects. Consistent results were 
obtained by RiceetaL using three-dimensional smoothed 

particle hydrodynamic (SPH) simulations, although the actual cool- 
ing ti me boundary was f ound to depend on the chosen equation of 
state (iRice. Lodato & Armitage 2005). 

Since self-gravitating discs transport angular momentum out- 
wards, it is useful to describe them as having an effective viscosity. 
IShakura & SunvaevI (Il973h suggested that an appropriate form for 
the viscosity in an accretion disc is 



aCsH, 



(3) 



where a « 1 i s the viscosity parameter and H is the disc 
thickness. iGammid (l200lh showed that the cooling time in a self- 
gravitating accretion disc could be related to the viscosity through 
an effective gravitational a that satisfies 



tteff = 



97(7 -1)/^' 



(4) 



where 7 is the s pecific heat ratio. It has since been shown 
dCossins. Lodato & Clarke 2009) that the amplitude of the pertur- 
bations in such a disc is related to aeff (or /3) through 

< > 1/2 1 

^E^""- oc— . (5) 

Essentially the current understanding is that the amplitude of the 
perturbations increase with increasing tteff and, hence, there is 
a maximum value for agff that a self-gravitating disc can sus- 
tain without f ragmenting. This maximum has been shown t o be 
Qma x ^ 0.06 taammiel |200ll : lRice. Lodato & Armitagll2005h . 

iMeru & Bate (201 1]) have, however, recently shown that simu- 
lations of self-gravitating discs using SPH fail to converge. In their 
simulations, the /3 value at which fragmentation occurs increases as 
the resolution increases. They do infer that this may be a numerical 
effect, but also suggest that this could mean that fragmentation may 
occur for arbitrarily long cooling times. If fragmentation requires 
short cooling times, then it becomes very unlikely that it ca n occur 
in the inner regions of protostellar discs CRafikov 2005 : Bolev et al.l 
I2OO6I : IStamatellos & Whitworth|[2008l : Iciarkell2009l) and, hence, is 
no longer a viable mechanism for the formation of gas giant plan- 
ets. On the other hand, if fragmentation can occur for very long 
cooling times it would mean that gas giants planets could form via 
gravitational collapse. However, this would be somewhat surprising 
as it would imply that fragmentation could occur for infinitesimally 
small perturbation amplitudes. We suggest, in this paper, that it is 
indeed a numerical effect and is related to the manner in which 
cooling is typically impl emented in SP H. 

That the results in iMeru & Bat3 (1201 ih are likely to be nu- 
merical suggests that our understanding of self-gravitating accre- 
tion discs is effectively unchanged. In a typical isolated disc that 
has a relatively low mass compared to the mass of the central star, 
the evolution is largely determined by Q and the cooling time. If 
Q ^ 1 and the cooling time is such that aeff < 0.06, it settles into 
a quasi- steady, marginally stable state with relatively small pertur- 
bations and in which angular momentum is transported outwards. 



If tteff > 0.06 the perturbations become non-linear and fragmenta- 
tion occurs. 

This paper is organised as follows. In Section 2 we describe 
SPH and how the implementation of the cooling may introduce nu- 
merical effects. In Section 3 we describe our results and show that 
implementing cooling in a manner that is more consistent with the 
SPH formalism may lead to convergence. In Section 4 we discuss 
these results and reach a few conclusions. 



2 SMOOTHED PARTICLE HYDRODYNAMICS 
2.1 The Basic Formalism 

Smoothed Particle Hydrodynamics (SPH) is a Lagrangian hydro- 
dynamic formalism in w hich a fluid is repres ented by pseudo- 
particles (see, for example jBenJ 199Ql : lMonaghanll 19921) . Each par- 
ticle in the simulation has a mass (m), position (x,y, z), velocity 
(Vx ,Vy,Vz) and internal energy per unit mass (u). 

Although SPH uses particles, the interpretation is that each 
particle represents a smeared out distribution of density. The contri- 
bution that particle j, located at r^, makes to the density at location 



(6) 



where W is the "smoothing kernel" that describes the form of the 
mass distribution and hj is the smoothing length associated with 
particle j. The density at r is then determined by summing over the 
N particles that contribute to the density at that location 



(7) 



To perform a basic hydrodynamic simulation, SPH needs an 
equation of state and has to solve the continuity equation, the mo- 
mentum equation and the energy equation. If the energy equation 
is included then typically an ideal gas equation of state will be used 



Pj = (7 - ^)puj 



(8) 



where 7 is the specific heat ratio, p is the fluid density at the lo- 
cation of particle j, and Uj is the internal energy per unit mass 
associated with particle j. 

The continuity equation is 

dp 



dt 



+ V • (pv) - 0, 



(9) 



but since SPH uses particles with constant masses, mass is con- 
served trivially and the continuity equation does not need to be 
solved explicitly. Ignoring self-gravity, the Lagrangian form of the 
momentum equation is 

5 = - -VP, (10) 

at p 

showing that - in the absence of gravity - the acceleration is deter- 
mined only by the local pressure gradient. In SPH, the acceleration 
of each particle is determined and Equation ([TO ]) ^i s typically rewrit - 
ten into a suitable form shown below ( Benz 199Ql : lMonaghanlll992h 



P^ 



(11) 



where the sum is over all the neighbours (i) of particle j and 
ViW{\ri — rj\;h) is the gradient of the kernel (a known func- 
tion) at the location of each neighbour, with h typically the average 
of the smoothing lengths for particles i and j. 
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The energy equation is 

de 

_ + V • (ev) = -PV • V, (12) 

where e = pu is the internal energy density. Using the continuity 
equation (Equation ^) the energy equation can be written in its 
Lagrangian form 

du _ P, 



-V-v 



(13) 



whic h can then also be written into a form more suitable for SPH 
re.g. JMonaghanl (Il992h ) 

^4EQ + |)v.-V.I^(|n-r,|;.), (14) 
where v^j is the velocity of particle i with respect to particle j. 



2.2 Introducing cooling 

To investigate how self-gra vitating discs evolve in the presen ce of 
cooling, many simulations f Gammi 3 l200ll :l Rice et al ] l2003 l) have 
used what is now referred to as /3-cooling and is illustrated in Equa- 
tion (|2]). The cooling term in the energy equation would therefore 
be 



dt 



e 

Tc 



(15) 



If this is introduced into the energy equation, Equation i\3l 
then becomes 

du P„ e 1 P. 



P^ e 1 
^ Vv 

dt p pre 



-V-v 

P Tc 



(16) 



Since each particle has an internal energy per unit mass (uj) asso- 
ciated with it, typically the cooling term for each particle is simply 
assumed to be 



duj 
~dt 



(17) 



This cooling form has been u s ed in many previous SPH sim- 
ulations (e.g. . Rice et al. 2003; Rice. Lodato & Armitagd [2OO5I : 
Cossins. Lodat o & Clarke 2009) and is the form used by 
Meru & Batel |2011) in their simulations that appear not to con- 
verge as the resolution is increased. 

The internal energy per unit mass associated with a particle is, 
however, the total amount of thermal energy that the particle has 
and, in an SPH sense - as with the mass - should be regarded as 
a smeared out distribution of thermal energy. We argue that using 
the standard /3-cooling formalism produces a mismatch in scales. 
The hydrodynamic equations are solved at the scale of the smooth- 
ing length h - if the gravitational potential is kernel- softened, then 
it can be argued that gravitational forces are also limited by the 
smoothing length h. For simulations that are not well-resolved, a 
density peak of leng th- scale A will be supp ressed by a factor of 
1 + /i/A, as noted bv lLodato & Clarkj (I2OI I), who use this to de- 
rive a resolution requirement for fragmentation: 



/3o{l + h/Hy 



(18) 



In implementing /3-cooling at the locations of particles only, 
the cooling operates far below the /i-scale. Purely hydrodynamic 
simulations will calculate in a smoothed sense, but adding an 
unsmoothed /3-cooling term violates this. We should then seek a 
means by which the /3-cooling term can be limited to the /i- scale. 
In SPH, any quantity. A, associated with the fluid can be estimated 



using densit y-weighted interpolations between the known particle 
values fe.g. jBodenheimer et 



A{r) = ^m.^Wi\r-r,\,h). 

.-1 Pi 



(19) 



Formally, the internal energy density at the location of particle j 
at position r - should then be determined using 

m 



e{r)=pir)y2—mW{\r-ri\,h). 



(20) 



The internal energy per unit mass at the location of particle j is 
therefore, formally, e(r)/y9(r). In principle, one would expect that 
- at the location of each particle - the interpolated value should be 
the same as the value assigned to the particle. This is the case in any 
regions that vary smoothly, but is n ot the case at shocks and other 
discontinuities ( Brownlee et aP fioOTil . That fragmentation tends to 
occur in the spiral shock regions may, therefore, suggest that using 
the unsmoothed thermal energies per unit mass may result in errors 
in the cooling rate that could be removed by using the smoothed 
thermal energy. We should acknowledge that there is a potential in- 
consistency in that the gas pressure is typically defined using Equa- 
tion ([8]) in which a smoothed quantity ( p) is combined with an un - 
smoothed quantity (u). As suggested bv lRitchie & ThomasI (I2OOII) . 
it may be more consistent to define the pressure as 



Pj = (7 - l)p(r) J2 — - r^l, h). 

i Pi 



(21) 



However, the pressure force and heating terms in Equations (fTTT) 
and ([T4I) are determined via interpolation and so it is not clear that 
the pressure itself has to be determined via interpolation. 

If the above is correct, then it is possible that the non- 
convergence seen in Meru & Bate! (l201lh is due to this mismatch 
between Uj and the smoothed internal energy per unit mass at the 
location of particle j. We could simply replace Equation fTTI ) with 



duj 
~dt 



I c Pt 



,h). 



(22) 



This, however, has a problem because the cooling is exponential 
and, hence, any particle cooler than all its neighbours will very 
quickly (and artifically) cool to a very low thermal energy. In- 
stead, we distribute the cooling, associated with particle j, across 
the neigbour sphere and use 



duj 
~dt 



1 W{0, hj)mjUj 



du 



dt 



1 rrii 



pj 



U^W{\YJ - Yi\,h), 



(23) 



cool 

where the lower of the two equations in Equation (l23l ) is applied to 
all neighbours i of particle j. Ultimately, the total cooling rate for 
particle j will be 



duj 
~dt 



lj2^u,W{\r,-r.lh), 
''c i Pi 



(24) 



which ensures that the correct amount of energy is removed from 
the neighbour sphere, but prevents the coolest particle from cooling 
exponentially to an artificially low thermal energy. We suggest here 
that using Equation (l24l) to represent /3-cooling is more appropriate 
than using Equation JTt] ) as it more correctly represents the cooling 
rate at the location of particle j and removes the correct amount of 
energy per time from the neighbour sphere, which is the smallest 
scale we can expect the SPH formalism to correctly accommodate. 
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2.3 Shock-tube test 

As disc ussed above, w e sug gest that the non-convergence seen 
in the iMeru & Bate! (l201lh simulations could be a conse- 
quence of their implementation of /3-cooling. As with many 
others (e.g., Riceetal. 2003, Rice, Lodato & Armitage 12005 



ICossins Lodato & Clarkell2009l) the cooling rate is determined us- 
ing the internal energy per unit mass assigned to each particle. We 
propose that it is more appropriate to use a smoothed thermal en- 
ergy (e.g.. Equation (|20l )) to determine the cooling rate. 

However, the interpolation used in SPH (e.g.. Equation ([19])) 
should, when applied at the location of a particle, recover the value 
of the fluid quantity assigned to that particle. If so, one would ex- 
pect that smoothed cooling and basic cooling should produce the 
same results. There are, however, regions - in particular at disconti- 
nuities such as shock waves or contact discontinuitie s - where this 
is not the case (see for example ferownlee et al ]l2007h . As a simple 
test we c onsidered a one-dimensional shock tube test problem (e.g., 
iHerquis t & Katz 1989) which we carried out usin g SPH and using 
the grid-based PENCIL CODE (lBrandenbur^l2003h . 

Figure [T] shows the density, p, (dashed line) and pressure 
(P = (7 — l)pu) divided by (7 — 1) (solid line) from the PEN- 
CIL CODE compared with the corresponding values from the SPH 
simulation (symbols). The pressure is divided by (7 — 1) simply 
to seperate it from the density. In the top panel it is clear that at 
the contact discontinuity (located at x ^ 0.1), the PENCIL CODE 
pressure is continuous (as expected) while the SPH pressure - de- 
termined using the unsmoothed internal energy per unit mass - is 
not. This discontinuity in the SPH pressure at the contact discon- 
tinuity is unphysical but, since the pressure force (Equation (HTI) ) 
and energy term (Equation ([14])) are determined by interpolation, 
the SPH code does not behave as if this discontinuity is actually 
present. If the cooling, however, uses the unsmoothed internal en- 
ergy, then this discontinuity will be present in the cooling and will 
produce an unphysical mismatch between cooling and heating near 
the contact discontinuity. Since fragmentation typically occurs in 
the region behind the shock, this mismatch could indeed play an 
important role in promoting fragmentation. The lower panel in Fig- 
ure [1] shows that the discontinuity in pressure at the contact dis- 
continuity is not present if the smoothed internal energy is used to 
determine the pressure. This suggests that using the smoothed inter- 
nal energy in the cooling formalism is more physically correct than 
using the unsmoothed internal energy. We therefore, below, carry 
out a series of simulations to determine if convergence is achieved 
when Equation (l24l) is used instead of Equation ( [TTl ). 



3 RESULTS 

AH of the sim ulations presented have the same bas ic conditions as 
iMeru & Bate! ([2OII1) , which are based on those of iLodato & Ricd 
(I2OO4I) . A central star with a mass of Mstar = 1 surrounded by 
a disc extending from — 0.25 to rout = 25, with a mass of 
Mdisc = 0. 1 and with a surface density profile of S oc . These 
are, as with lMeru & Batd (1201 ih . scale-free simulations and so we 
don't explicitly assume any scale, although typically these could 
be regarded as representing a disc from 0.25 to 25 au around a star 
with a mass of 1 M© . 
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Figure 1. Comparison between an SPH shock tube test problem (symbols) 
and one performed using the PENCIL CODE (lines). These show a shock 
wave at x ~ 0.3 and a contact discontinuity at x ~ 0.1. In the top panel, 
the SPH pressure (open circles) is determined using the unsmoothed in- 
ternal energy resulting in an unphysical jump at the contact discontinuity. 
Using the smoothed internal energy (lower panel) removes this unphysical 
discontinuity in the pressure. 



3.1 Comparison with earlier simulations 

To ensure that our SPH code is comparable to that of lMeru & Batd 
(>2011) we first run simulations using the basic cooling as described 
by Equation ( [TTl ). Table[T] shows the simulati ons that we performe d 
using basic cooling. You'll notice that, unlike Meru & Bate" ('2011'), 
we do not carry out an y 31250 particle simulations. As shown by 
iBate & Burkerj (Il997h . an SPH simulation only correctly repre- 
sents fragmentation if the minimum resolvable mass 



2Mu 



neigh 



iVtc 



(25) 



is less than the Jeans mass. In Equation ([25l) Mtot is the total disc 
mass, A^neigh is the number of neighbours each particle has (typi- 
cally 50) and Atot is the total number of particles in the simulation. 

The Jeans mass can be approximated using Mj — ttT^H^ 
where H = Cs/^ is the disc scale height. If we use H = ^o/r 
then, given that Mdisc = 0.1, So = 0.1/(27r[rout - Hn]) ^ 
0.1/(27rrout). Since th ese disc s are self-gravitating, we can use 
Q CsVL/{iiGi:) 1 (lToomra.1964) to give 



1/2 



(26) 
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No. of portloles 

Figure 2. Figure showing, for an SPH simulation with E oc r~^, Mdisk = 
0.1, Mstar = 1, and rout = 25, the radial range where the Jeans mass is 
resolved (unfilled) and unresolved (filled) plotted against particle number. 
This clearly shows that if fewer than 80000 particles are used, the Jeans 
mass is not properly resolved anywhere in the disc. 



Table 1. List of the simulations using basic cooling. 



Simulation No. of particles ^ Fragment? 



250000 
250000 
250000 
250000 
500000 
500000 
500000 
2000000 



Yes 
Yes 
No 
No 
Yes 
Yes 
No 
Yes 



where we've also used that Mstar = 1 and G — 1. We can therefore 
show that the Jeans mass in these simulations has the following 
dependence on radius 



Mj 



O.OOlr'' 



(27) 



Figure [2] shows a plot of radius against particle number where the 
filled region shows the radial range for which the Jeans mass is not 
resolved while the unfilled region shows the radial range for which 
the Jeans mass is resolved. Figure [2] shows clearly that if fewer 
than 80000 particles are used, the Jeans mass is not resolved any- 
where in the disc. With 250000 particles, the Jeans mass is resolved 
beyond r — 11 and for 2 million pa rticles it is resolved beyond 
r = 8. That the iMeru & Bat3 (1201 ll) results include simulations 
with 31250 particle, and that these simulations are consistent with 
the resulting trend, suggests that numerical effects are influencing 
the results. We also don't carry out any simulations with more than 
2 million particles simply because our code is parallelised using 
OpenMP and we are therefore limited as to how many processors 
we can access and, hence, can't realistically carry out higher reso- 
lution simulations. 

Table [T] shows tha t these simulations are consistent with those 
of lMeru&Bat9 ( l201lh . For 250000 particles, fragmentation occurs 
for 13 between 5 and 6, for 500000 particles it is between 7 and 8 
and for 2 million particles fragmentation occurs for /3 > 8. Figs. 
[SJHI and[5]show surface density images of a number of the simula- 




Figure 3. Simulations with 250000 particles and using basic cooling. The 
left panel is /3 = 4 while the right panel is /3 = 6. 




Figure 4. Simulations with 500000 particles and using basic cooling. Left 
panel is /5 = 7 and the right panel is /3 = 8. 



tions. For the 250000 and 500000 particle simulations we show an 
image of simulation that fragmented and one that did not. For the 
2 million particle simulation we simply illustrate (in Fig. [5]) that 
the /3 = 8 simulation did indeed undergo fragmentation. These 
simulations are all run for at least 6 outer rotation periods and, if 
there is any evidence for fragmentation, they are continued until 
either these fragments become much denser than the local density, 
or they shear away. The clumps in Figs. [3] |4] and [5] are therefore 
smaller and fewer in number than would be the case if the image 
showed an epoch just after fragmentation started. 



3.2 Smoothed cooling 

Having shown that, using the basic co oling represented by Equation 
(ITtI ), we get results consistent with Meru & Bate! (|201 ih we then 
carried out a large number of simulations using what we shall refer 




Figure 5. Simulation with 2 million particles, using basic cooling and with 
/3 = 8. 
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Table 2. List of the simulations using smoothed coohng. 



Simulation 


No. of particles 




Fragment? 


\ 


950000 


4 


Yes 


2 


950000 


4.5 


Yes 


3 




5 


No 


4 


250000 


6 


No 


5 


250000 


7 


No 


6 


500000 


5 


Yes 


7 


500000 


6 


Yes 


8 


500000 


7 


No 


9 


500000 


8 


No 


10 


2000000 


5 


Yes 


11 


2000000 


6 


Yes 


12 


2000000 


7 


No 


13 


2000000 


8 


No 


14 


10000000* 


7 


Yes 


15 


10000000* 


8 


Yes 


16 


10000000* 


9 


No 




Figure 6. Simulations with 250000 particles and using smoothed cooling. 
The left panel is /5 = 4, and the right panel is /3 = 5. 



to as smoothed cooling (Equation (|24])). Table |2] lists all of these 
simulations and whether or not fragmentation occured. The starred 
10 million particle simulations are ones in which 4 million particles 
(with a total mass of Mdisc = 0.04) were placed in a ring from 
r = 15tor = 25. This is equivalent to 10 million particles (with a 
total mass of Mdisc = 0.1) placed between r — 1 and r = 25, and 
so we refer to these as 10 million particle simulations. 

Figs. [6l |7] [8] and |9] show surface density images from the 
simulations using smoothed cooling showing the range of /3 val- 
ues across which the simulations go from fragmenting to non- 
fragmenting. 

The results from iMeru & Bate! (l201l l) suggest that the de- 
pendence of the fragmentation boundary (/3crit) on resolution is 
/3crit oc n^/^ with n the number of particles in the simulation. 
This is equivalent to saying that the fragmention boundary depends 
linearly on resolution (i.e., h oc n^^^). Table [2 immediately shows 
that ou r results are considerably different to those of Meru & Bate 
(I20T ll) with the fragmentation boundary varying from between 4.5 
and 5 for simulations with 250000 particles to between 8 and 9 for 
the simulation with effectively 1 million particles. Fig. [Tol com- 
pares the results from'Meru & Bate' ('2OII') with the results obtained 
here. The filled symbols are from Meru & Bate ( 201 1) and the open 
symbols are from this work. Triangles represents simulations that 
fragmented whi le squares are for tho se that did not. The filled cir- 
cle represents a lMeru & Bate! (l201lh simulation that they referred 
to as being borderline. We see very little evidence for what might 




Figure 7. Simulations with 500000 particles and using smoothed cooling. 
The (5 values are /3 = 4 (top left), 5 (top right), 6 (bottom left), and 7 
(bottom right). 



Figure 8. Simulations with 2 million particles and using smoothed cooling. 
The top4eft panel is ^ = 5, the top right-panel is /3 = 6 and the bottom 
panel is /3 = 7. 



be regarded as borderline systems. Our simulations either definitely 
fragment (the clumps continue to grow and become extremely 
dense) or the system settles into a quasi- steady state with no signs 
of fragmentation. In some cases clumps will form and shear away, 
but this appear s to happen v ery rarely. We would generally inter- 
pret the lMeru&Batd (1201 ih borderline cases as non-fragmenting 
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Figure 9. Simulations with a ring consisting of 4 million particles and using 
smoothed cooling. This is equivalent to a disc simulation with 10 million 
particles. Left panel is /3 = 8 and the right panel is /3 = 9. 
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Figure 10. Figure showing fragmentation boundary from Meru & Bate 
(filled symbols) together with results from this work (open symbols). The 
line repre sent the relat i onship between /3crit and particle number deter- 
mined bv Meru & Batd 120111) . 



systems. The line in Fig. [TOl show s the relationship bet ween /3crit 
and particle number determined bv'Me ru & Batd (1201 ih . 

It is clear from Fig.[TO]that using smoothed cooling (Equation 
(|24] )) produces results significantly different to that produced by 
simulations using basic cooling (Equation (fTvl)). Although the ef- 
fective 10 million particle simulation has a fragmentation boundary 
slightly higher than the 2 million and 500000 particle simulations, 
it does appear as though we can conclude that the fragmentation 
boundary may be converging to a value of /3crit < 10. That the 
500000 particle and 2 million particles simulations produced very 
similar values for /3crit (between [3 = 6 and [3 = 7) had lead us 
to initially conclude that these simulations had already converged. 
The higher value for /3crit produced by the 10 million particle sim- 
ulation does suggest that convergence hasn't yet been reached, but 
could also imply that a ring of particles is not an exact proxy for a 
more massive and larger disc with more particles. 

3.3 Comparison of cooling rates 

That the interpolated value of the thermal energy should approxi- 
mate well the individual particle values everywhere except at dis- 
continuities suggests that there shouldn't be a significant differ- 
ence, in general, between the cooling rate in simulations using ba- 
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Figure 11. Comparison of the cooling rate using smoothed cooling (solid 
line) and basic cooling (dashed line). These cooling rates are azimuthally 
averaged and also averaged over the final 1.5 outer rotation periods. The 
inner regions of the disc is dominated by numerical viscosity which artifi- 
cially heats this inner region. Beyond r = 5, there can be local variations 
that are as large as 10 %, but globally the two cooling formalisms are very 
similar, differing by only a few percent. 



sic cooling when compared to those using smoothed cooling. To 
test this we considered azimuthally averaged cooling rates for 2 of 
the simulations with the same number of particles (500000) and 
both with [3 = 8 that did not fragment when either smoothed or 
basic cooling was used. We also averaged over the final 10 dumps 
covering a time of 1.5 outer rotation periods. The result is shown in 
Figure [TT] with the solid line being the cooling rate using smoothed 
cooling and the dashed line being the cooling rate using basic cool- 
ing. 

Figure [TT] show that the cooling rates are indeed similar and 
have the same radial trend. The large increase in cooling rate in- 
side r = 5 is a consequence of this region being poorly re- 
solved due to particles having been accreted onto the central star 
(Bate, Bonnell & Price 1995), resulting in the numerical viscos- 
ity dominating in this inner region (see for example Forgan et all 
2011). This artificially heats the inner disc and produces the large 
increase in cooling rate. Beyond r = 5, where numerical viscos- 
ity does not dominate and the system is well resolved, the local 
variations can, however, be as large as 10 %. If, however, we sum 
the individual particle cooling rates from a single dump and take 
the average for r > 5, the smoothed and basic cooling formalisms 
differ by, typically, less than 5 % with neither cooling formalism 
being, consistently, the most efficient. If this is averaged across the 
final 1.5 outer rotation periods (10 dumps) the difference is only 1.8 
%. This suggests that even though locally there can be reasonable 
(10 %) differences between the two cooling formalisms, globally 
they are very similar. 



4 DISCUSSION AND CONCLUSIONS 

To study the evolution of self-gravitating accretion discs that are 
undergoing cooling, many authors have used a form of cooling 
known as /3-cooling and which is represented by Equation (|2]). In 
SPH, this has generally been implemented by cooling each particle 
with a cooling rate determined using the internal energy per unit 
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mass associated with that particle. iMem & Bate! (1201 ih have re- 
cently shown, however, that such simulations do not converge and 
that the value for /3 at which fragmentation occurs (/3crit) increases 
linearly with resolution (i.e., /3crit oc n^^^). 

We confirm here that this non-convergence is indeed seen in 
such simulations, but suggest that it may be a SPH problem and 
may result from using the internal energy per unit mass associated 
with each particle. This presents a mismatch of scales between the 
/i- scale (where the traditional SPH equations operate) and the es- 
sentially infinitesimal scales at which the additional /3-cooling term 
operates. We advocate an alternative approach to /3-cooling where 
the cooling rate is determined by smoothing over contributions due 
to the neighbours. This correctly represents the local internal en- 
ergy per unit mass, and removes any unphysical discontinuities (in 
pressure for example) that could occur in the shocked regions when 
using the unsmoothed internal energy. Although our simulations 
using this smoothed cooling do not strictly converge, they appear 
to be converging to /3crit < 10 (for 7 = 5/3), consistent with re- 
sults from earlier work. To further test the convergence of our pro- 
posed cooling method would require higher resolution simulations 
(> 10 million particles) that we are not able, at this stage, to carry 
out in a reasonable time. We also compare the cooling rates using 
these two cooling formalisms. Although there are local variations 
that can be as high as 10 %, globally (if we consider the regions be- 
yond r = 5, that are well resolved and not dominated by numerical 
viscosity) the two cooling methods differ by only a few percent. 
This is consistent with our suggestion that it is local differences in 
the two cooling methods (occuring primarily at shocks and other 
discontinuities) that result in different fragmentation boundaries. 

^It has also been suggested (j Paardekoo per, Baru teau & Merul 

l201lh that this non-convergence is due to the use of smooth initial 
conditions. Fragmentation then occurs at the boundary between the 
turbulent and laminar regions as the system cools towards a quasi- 
steady state. It is suggested that it may be more appropriate to pro- 
duce an initial quasi- steady state using a long cooling time (a high 
13 value) and then to vary the value of (3. This is a very reasonable 
suggestion but comparisons of their work with this paper may not 
be wholly appropriate. Their implementation of /3-cooling in a 2D 
grid is closer to our smoothed cooling techniques than traditional 
SPH /3-cooling, which (given the results of this paper) would assist 
efforts to converge. On the other hand, their simulations require a 
seed of white noise to be added to the smooth initial conditions to 
generate spiral structure. SPH simulations typically contain Pois- 
son noise in the initial conditions as a result of innate particle dis- 
order, and the resulting fragmentatio n is rel atively insensitive to the 
level of this noise (Cartwright etal ]|200i). These differing noise 
properties may provide a reason as to why grid-based simulations 
converge more easily after being relaxed. In any case, we argue that 
their results are n ot sufficient to explain the linear dependence on 
resolution seen bv lMeru&Batd ( l201ir) . 

Although we have focused on simulations that use {3- 
cooling, the results presented here also have potential implica- 
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20091) . In most of these formalisms, the tem- 



perature at the location of a particle is determined from its internal 
energy per unit mass. It is possible that it too should be determined 
by smoothing across the neighbour sphere. The form of the cooling 
in these radiation transform formalisms is, however, very different 
to that for /3-cooling and typically these simulations cool towards 
an equilibrium temperature rather than exponentially towards zero 



(as is the case for /3-cooling), so the convergence problem seen with 
/3-cooling may not affect these simulations. 

Ult imately what w e suggest here is that the non-convergence 
seen bv iMeru & Bate! ([2OII) is primarily a numerical effect and 
hence that there is no evidence to suggest that fragmentation can 
occur for arbitrarily long cooling times in self-gravitating accretion 
discs. The standard picture that fragmentation requires short cool- 
ing times (tteff > 0.06) still, therefore, appears to be valid. This 
is c onsistent \yith the results of Cossins, Lodat o & Clarkd (l2009l) 
and iRice et all (l201lh who show that the perturbation amplitudes 
in self-gravitating discs depend inversely on the cooling time and 
hence fragmentation occurs when the perturbation are large (non- 
linear) and, consequently, when the cooling time is short. 
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